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The viscous, Navier-Stokes solver for turbomachinery applications, 
MSUTC has been modified to include the rotating frame formulation. The three- 
dimensional thin-layer Navier-Stokes equations have been cast in a rotating 
Cartesian frame enabling the freezing of grid motion. This also allows the flow- 
field associated with an isolated rotor to be viewed as a steady-state problem. 
Consequently, local time stepping can be used to accelerate convergence. The 
formulation is validated by running NASA’s Rotor 67 as the test case. Results 
are compared between the rotating frame code and the absolute frame code. The 
use of the rotating frame approach greatly enhances the performance of the code 
with respect to savings in computing time, without degradation of the solution. 
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CHAPTER I 


INTRODUCTION 

The flow field associated with turbomachinery applications is very com- 
plex. This flow field is highly unsteady, involving a wide range of time and 
length scales [1]. The goal of a CFD code for simulating such flows is to be able 
to model all aspects of the flow phenomena reasonably well. The computing time 
that is required should also be practically feasible for the CFD code to be useful 
as a design tool for turbomachinery applications. 

The MSUTC code, which is a turbomachinery flow analysis research tool, 
is an attempt to achieve this goal. This code has been under development at the 
NSF ERC for Computational Field Simulation at Mississippi State University. 
It is a viscous flow solver which is capable of simulating even or uneven-blade- 
count, single-rotating, counter-rotating or rotor-stator, axisymmetric or non- 
symmetric, multistage, at angle-of attack geometries. Chen [2] noted that 
though converged unsteady solutions can be obtained using the above code, they 
are very expensive, because the use of Newton sub-iterations requires tremen- 
dous computing time. Thus, the performance of the code can be enhanced if one 
can avoid the Newton sub-iterations. 

Navier-Stokes equations which are the governing equations of fluid me- 
chanics can be expressed in the vector invariant form so that they are indepen- 
dent of the coordinate system. For flows in rotating machinery, it is convenient 
to cast the equations in the rotating coordinate system. Adamcyzk et al. [3] cast 
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the equations in a rotating cylindrical coordinate system to simulate viscous 
flow through turbines. Although unsteady computations are very important in 
turbomachinery simulation, many of the problems can be viewed as steady-state 
problems in the rotating frame, e.g. an isolated rotor. Consequently, local time 
stepping and multigrid methods can be utilized to accelerate convergence. 

Due to this motivation, at the start of this study an effort was made to de- 
velop a cylindrical coordinate system Navier-Stokes code for simulating the 
steady flow through a single blade row machine. During the course of the inves- 
tigation it was realized that only a few modifications were necessary to incorpo- 
rate the rotating frame formulation into the existing unsteady code if one re- 
tained the Cartesian coordinate system but solved the problem in a rotating 
frame similar to the approach in [4] and [5]. This is the approach that was 
adopted for the work presented in this study. The existing viscous solver is modi- 
fied to include the rotating frame formulation using the absolute velocity vector. 

This thesis is divided into five chapters. Chapter II develops the govern- 
ing equations cast in the rotating Cartesian coordinate frame. The velocity used 
in the equations, however is the absolute velocity represented in the rotating 
frame. The equations are then non-dimensionalized and transformed from Car- 
tesian coordinates to steady curvilinear coordinates to enable the use of a body 
conforming grid. The thin layer approximation is also explained in this chapter. 
Chapter III describes the aspects of the flow solver with all the changes that are 
made due to the change in the formulation. The implicit formulation, Newton’s 
method, flux formulation, flux Jacobians, the N-Pass solution scheme and the 
boundary conditions are all explained and discussed in this chapter. Chapter IV 
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presents the results obtained from the modified code. The test case used is 
NASA’s Rotor 67. The results obtained from the rotating frame approach are 
compared with those obtained form the absolute frame approach. Chapter V is 
devoted to summary and conclusions. 


CHAPTER II 


GOVERNING EQUATIONS 

The development of the Navier-Stokes equations in the rotating coordi- 
nates are presented in this chapter. The primary objective in this study is to 
make use of the quantities in the fixed frame formulation so that a minimum 
amount of changes need to be made to the existing code. To achieve this goal, the 
formulation presented in this chapter develops the Navier-Stokes equations in 
rotating coordinates using the absolute velocity, i.e. velocity with respect to the 
fixed inertial frame. 

The differential forms of the equations of conservation of mass, momen- 
tum and energy can be collectively stated in a general form. The differential 
form of the general conservation principle can be stated as [6] : 

# + div f = C (2.1) 

dr 

where f = Au + B and A, B, C, are tensor quantities such that A and C 

have the same tensorial order, and if B 0 then it is an order higher than 

A and C. Equation (2.1) is written with respect to an inertial frame so that the 
vector u is the absolute velocity. Consider a transformation from the inertial 

frame to a moving non-inertial frame which is of the form [6]: 
x l = x l (x 1 ,x 2 ,x 3 ,t), i = 1,2,3 and r = t 
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where (x v x 2 , x 3 ) are the rectangular Cartesian coordinates defining the some- 
how, inertial frame, (at 1 , x 2 , x 3 ) are the general coordinates defining the mov- 
ing frame and a are the covariant base vectors of the moving frame. If the rela- 

tive velocity of a fluid particle, which is the velocity with respect to the 
nonsteady coordinates, is denoted by v and the velocity of the moving coordi- 
nates is denoted by w ; the absolute velocity of a fluid particle can be expressed 
by the following equation [6] : 

u = v - w (2.2) 


where 


w 



(2.3) 


The contravariant components of w are the partial time derivatives of each of 

the coordinates in the moving frame and it can be proven to be negative of the 
velocity of the moving frame with respect to the inertial frame. Thus, the vector 
w can be thought of as the velocity of the absolute frame, which an observer in 

the rotating frame would notice. Details of the derivation are presented in Ap- 
pendix A. 

Due to the above transformation, the following relation exists between the par- 
tial time derivatives in the two frames of reference [6]: 


dj ) 

dt 


d( ) 
dr 


-I- w • grad( ) 


(2.4) 


Thus, the unified conservation law, equation (2.1), transforms to [6]: 
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+ (grad A) • w + div f = C (2.5) 

or _ 

The rotating Cartesian coordinate frame utilized in this study is illustrated be- 
low. 



z 


Figure 2.1 Rotating Cartesian Coordinate Frame 

It is evident from Figure 2.1 that the rotating frame is rectangular Carte- 
sian in nature, due to which there is no distinction between covariant and con- 
travariant components. Also the base vectors do not vary in space, so that: 

da 

- 4=0 ( 2 . 6 ) 

dx l 

The temporal derivatives of the base vectors are given by [6]: 
da dw 

^1+^=0 (2.7) 

dr dx l 

Choosing a rotating frame illustrated in Figure 2.1 results in: w = - Q x r 

where Q - - Qa . Therefore, w = - Qx 3 a + Qx 2 a . Consequently, 

_i _ -2 ~3 
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dw 

~dx l 


j- _ i dwJ 

div w = — - • a = 


a 


dx l ~j 

Using the following tensor identity : 


a i _ dwi. _ o 


dx l 


(2.8) 


div (Aw) = (grad A) • w + (div w)A = (grad A) • w 


(2.9) 


and equation (2.8), one can rewrite equation (2.5) as: 

^ + div (Au;) + div f = C (2.10) 

As explained at the beginning of this chapter, the primary objective is to 
cast the Navier-Stokes equation in the rotating frame using the absolute veloc- 
ity vector. The continuity, momentum and energy equations that are derived be- 
low in this chapter, have used the absolute velocity vector u. However, all the 

vectors appearing in equation (2.10) are represented using the base vectors of 
the rotating frame. This concept can appear confusing. This confusion is re- 
moved by addressing two important aspects related to any velocity vector in the 
above equations. The first being definition of the vector, and the second is repre- 
sentation of the vector. 

The keywords to defining a velocity vector are: with respect to. This clari- 
fies which reference frame one is referring to when defining the velocity vector. 
Thus, the absolute velocity u, is defined as the velocity of a fluid particle in the 

concerned flowfield, with respect to the fixed inertial frame. Therefore, a station- 
ary observer who selects a right handed rectangular coordinate system like the 
one in Figure 2.1, but one which is not moving, observes the velocity of a fluid 
particle in the flow field as u. On the other hand, an observer moving with the 
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axes described in Figure 2.1, observes the velocity of the fluid particle in the flow 
field as v . Thus, the relative velocity v , is defined as the velocity of a fluid par- 
ticle with respect to the rotating non-inertial reference frame. 

Once the definition of the vector is taken care of, the issue of representing 
the vector arises. These two issues are independent of each other and once a vec- 
tor is defined with respect to a frame, there is no need to express the components 
of the same vector using a coordinate system in the same frame. Thus, the abso- 
lute velocity vector can be expressed in component form with respect to the basis 
of the rotating Cartesian coordinate system shown in Figure 2.1. Since the exer- 
cise in this chapter is to cast the Navier-Stokes equations in the rotating frame, 
all the mathematical operations arising in the equations below, viz., time deriv- 
atives, curl, grad, div, are carried out in the rotating frame and involve its base 
vectors. To be consistent and to simplify mathematical operations, the rotating 
frame has been chosen to represent all vectors and tensors. 

2.1 Continuity Equation 

The continuity equation is obtained from the generalized conservation 
law, equation (2.10), by substituting A = Q, f = qu, B = 0, C = 0. Thus, 
the continuity equation in the rotating frame can be written as: 

+ div (pw) + div (qu) = 0 (2.11) 

2.2 Momentum Equation 

The momentum equation in the non-inertial frame is obtained by substi- 
tuting A = qu, B = - T = pi - a, C = 0in equation (2.10). 
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The resulting equation is: 
b(gu ) 

— + div ( guw ) + div ( guu ) = div ( T ) (2.12) 

oT ~ _ 

where T is the stress tensor, p is the thermodynamic static pressure, I is the 

identity tensor and a is the deviatoric part of the stress tensor containing the 
viscous shear stress. It can be represented by: 

o = A(div u)I + //[(grad u ) + (grad i/) 7 ! (2.13) 

where// and A are the first and second coefficients of viscosity, respectively. Writ- 
ing u = u l a and utilizing the result in equation (2.7), the momentum equation 

~ — i 

can be rewritten as: 

a/ s\ ~ dw 

— f — a + div (guw) + div (quu) = div (T) + gu l — =r (2.14) 

oT — — dx l 


2.3 Energy Equation 

The generalized conservation law, equation (2.10), yields the energy 
equation by substituting A = e t , B = — T ■ u + q, C = 0, where 


e t = -At + ^g | u | 2 is the total energy and q is the conductive heat flux vector. 
y — X £ 


The energy equation in the rotating frame, therefore, is: 

be, 


dr 


7 + div (e t w) + div (e t u) = div (T ■ u) — div(q) 


(2.15) 


2.4 Conservation Law Vector form of the Equations 


Equations (2.11), (2.14) and (2.15) can be expanded and written in the 
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conservation law form by carrying out the mathematical operations of diver- 
gence and dot products on the various vectors and tensors involved. While doing 
the operations involving the base vectors one utilizes the results that the base 
vectors do not vary in space and that since the coordinate sixes are rectangular 

Cartesian in nature; a ■ a = d t . For simplicity, let 

~i ~j 

x 1 = x, x 2 = y, x 3 = z, x = t and u 1 = u, u 2 = v, u 3 = w. The five equa- 
tions that result can then be combined into one single vector equation of the 
form: 


df dg Qh df 3 ■ d g v . dh v . „ 
dt dx dy dz dx dy dz 


( 2 . 16 ) 


where the dependent variable vector q, the flux vectors f, g, h, the viscous flux 
vectors f°, h v , g v , and the source term vector s are defined as follows: 


<7 = 


' Q ' 

QU 

QV 

QW 


( 2 . 17 ) 


/• = 


QU 


QV ~ qQz 


QW + Q@y 

QU 2 + p 


qvu — quQz 


qwu + QuQy 

QUV 

, g = 

qv 2 + p — qv£2z 

, h = 

qwv + QvQy 

QUID 


qvw - qwQz 


qw 2 + p + QwQy 

u(e t + p) 


v(e t + p) - efiz 


w{e t + p) + efiy 

r 

0 

1 



f° = 


^XX 

r xy 

?xz 


UTxx UT'xy UJTxz Qx 


( 2 . 18 ) 
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g v = 


UT 


yx 


0 

Tyx 

T yy 

tyz 

+ VXyy + WTy Z + q y 


h v = 


U T, 


0 

Tzx 

r zy 

X Z z 

+ vr zy + wr zz 


+ Qz 


s = 


‘ 0 ' 
0 

— qwQ 
qvQ 

0 


(2.19) 


(2.20) 


The viscous stress terms appearing in the equation (2.19) are defined below un- 

2 

der the assumption that bulk viscosity is negligible, so that A = — . 


- f" (2 i - 

T = 2 (2 du_ 

r yy sy 

- §"< 2 f - 


dv 

dy 

du 

dx 

du 

dx 


dw ) 
dz ’ 

dW \ 

dz h 

dv\ 

dy’’ 


r _ _ _ ,,/dU , dVs 

*xy - r yx - fi( dy + dx ) 

T _ T _ ,, ( dv , dWs 

Tyz - r zy - M az + dy ) 

_ _ _ ../du; , dux 

t Z x - r xz - fi( dx + dz ) 


( 2 . 21 ) 


The conductive heat flux terms which appear in the viscous flux vectors are de- 
fined below. The conductive heat flux vector is defined as: q = — fc(grad T) 

where k is the thermal conductivity, and T is the temperature. The components 
of the heat flux vectors along each of the coordinate directions, therefore, are: 
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Qx 

Qy 

Qz 





( 2 . 22 ) 


All the components of the tensors and vectors that appeared in the equa- 
tion (2.16) are in terms of the base vectors of the rotating frame. The existing 
code uses a fixed frame formulation so that the components involved are in 
terms of the base vectors of the fixed frame. For the purpose of validation of the 
rotating frame formulation one needs to be able to compare the flow fields gener- 
ated as the solutions from the two formulations. Thus, a relation between com- 
ponents in the rotating frame and components in the fixed frame needs to be ob- 
tained. Letting u a , v a , and w a be the components of the absolute velocity with 
respect to a Cartesian coordinate system in the fixed frame, and, assuming for 
convenience, that the initial configuration of the rotating frame matches with 
the stationary frame (say at instant t = 0); the following relations are obtained: 
u a = u 

v a = vcos(Qt) - wsin(Qt) (2.23) 

w a = u sin {Qt) + wcos(Qt) 


where, it is emphasized that u, v, and w are the components of the absolute 
velocity with respect to a rotating Cartesian basis. It is apparent from the above 
equation that at instant t = 0 and after every subsequent full rotation, the com- 
ponents of the absolute velocity in the rotating frame are identical to the compo- 
nents in the stationary frame. 

Now, the Navier-Stokes equations have already been cast in the moving 
frame, so that there is no need to move the grid anymore. The flow field that is 
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generated in terms of the components, by solving the Navier-Stokes equations 
in the form presented in equation (2.16), for purposes of comparison, can be 
changed according to equation (2.23) to match the flow field which results from 
solving Navier-Stokes equations derived in the inertial frame. However, the dis- 
tinct advantage of the former over the latter is that local time stepping can be 
used for steady flows in the rotating frame. Another advantage is the freezing 
of grid motion because in the fixed frame formulation the grid motion is impera- 
tive. 


2.5 Non-Dimensionalization of the Equations 


All the variables used in the formulation till now are dimensional vari- 
ables. In order to obtain a form of the equation (2.16) with non-dimensional 
quantities, a scaling of all the dimensional variables (denoted by ~ ) is carried 
out using the following relations: 


X = 4, y = TT, 2 = 4, U = v=-r~, w= ] r-, Q = 

/) n n a n On a n 


DQ 


D 


D 


p _ 6 a o* 0 _ _e_ h = Jl „ = f* 

P — 2’ C? ~ A > £ A > e A A 2 ’ ™ A 2 9 ^ 

p 0 ao D Qo a o 


,2 ' 

ao 




T y a a 2 ’ T wall 
Q o a o 


= -)fTT. 9,, =- L£ 3, T = 4-, Zo = yRT 0 


^ogo 

D 


ix 

A A 


Q o°o 


The scaling has been accomplished using the reference quantities (de- 

A 

noted by subscript 0) which are the total conditions and the reference length D 
which is the maximum diameter of the blade tip. The detailed development of 
the scaling of the equations can be found in [2]. The scaling results in a form of 
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the vector equation which is exactly similar to equation (2.16) but with non-di- 
mensional variables. The form of the dependent variable, flux, viscous flux and 
the source term vectors is also exactly same as shown in equations (2.17) - (2.20), 
but with all variables non-dimensional. The only difference arising out of the 
scaling of equations is the introduction of two new non-dimensional quantities, 
viz., Reynolds and Prandtl numbers (Re, Pr), which appear in the non-dimen- 
sional form of the viscous shear stress and the conductive heat flux terms. The 
non-dimensional viscous shear stress and the conductive heat flux terms are: 


_ _ 1 2.., n du dv dUK 

Txx Re3^ dx by dz’' 


1 2 , 


dv dll dWs 


— ± — ulO _ ULU \ 

yy Re 3 by dx dz’' 

r - JL2../o dw du dV\ 

zz Re 3 ^ dz dx by’' 


_ - -r - 1 „( dU J. 

r X y - Ty X - + zJ 


dx' 


dv . dw % 


-r — -r — X tf( VU 4- UUJ \ 

Tyz ~ Tzy ~ + Q y > 

_ _ _ _ 1 „ ( dw , du-, 

Tzx - r xz - + dz ) 


(2.24) 


0 - J_ 


dT 

qx ~ Re 

_(y - l)Pr_ 

dx 

a - -L 


dT 

Qy ~ Re 

(y - l)Pr 

dy 

a — _ 1 _ 


dT 

qz Re 

(y - l)Pr 

dz 


(2.25) 


where Re = 


fio 


and Pr = 


c p fl 

A 

k 


2.6 Curvilinear Coordinate Transformation 

In order to carry out numerical computation of the governing equations 
on a configuration comprising of a complex geometry, it is simpler to use a body- 
conforming grid in which the solid surfaces in the geometry correspond to 
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constant coordinate surfaces. To this end, a curvilinear coordinate transforma- 
tion needs to be introduced in the rectangular Cartesian space being used in the 
formulation thus far. The following general, nonorthogonal, steady, curvilinear 
coordinate system is introduced in the ( x , y, z) space: 

£ = £(z>y*z) 

rj = r}(x,y,z) (226) 

C = t(x,y,z) 

T = t 


Only the final results due to the transformation are stated in this section. A de- 
tailed presentation of all the steps involved and the how the relations for the 
metric terms are obtained, is made in Appendix B. Using chain rule yields: 


d 

= d_ 




dt 

dr 




_d_ 

dx 

= t A 

+ 

^ x drj 

+ £ — 
Qx d$ 

_d_ 

dy 

Ms> 

ii 

+ 


+ t: — 

d_ 

dz 

= t A 
§2 d£ 

+ 

•>4 

+ £ — 


(2.27) 


The Jacobian J of the inverse transformation is: 
\d(x,y,z)\ 


J = det 


d(£,»7,C)l 

= x %(y>? z !; ~ z vy^ — y tf-X-qZ^ ~~ z t/X^) + z ^x^ — y rjX^) 


(2.28) 


The metric terms in equation (2.27) are: 

lx = - Zrjy c ) Vx = ~ Zx = - Z^rj) 

£>y — ~j^Z r pC XjjZf-) T]y — ZgpCf?) £_y J^z^rj (2.29) 
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Substituting relations obtained in equation (2.27) into equation (2. 16) re- 
sults in the form of the governing equations in stationary curvilinear coordi- 
nates. Only the final equation is stated here. The resulting form is: 


dQ.dF.dG. dH dF v 
dr dr] 


dr, ac 


+ S 


where 


m 


Q = J 


QU 

QV 

gw 

et 


K = J 


’ QK' ' 
guK' + k x p 
gvK' + k y p 
gwK' -I- k z p 
eJT + pK 


where K' = k x u + k y u + k z w + k t 
k t — — k y Qz + k z £2y 
K = k x U + kyV + k z w 


(2.30) 


(2.31) 


(2.32) 


r = j 


0 ' 

T kx 

T ky 

Tkz 

Q k 


s = J 


o ‘ 
0 

- gwQ 
gvQ 
0 


(2.33) 


(2.34) 


In the above equations (2.32) and (2.33), K = F, G, H\ K' — U’, V', W'; 
K = U, V, W; It = F°, G v , H v ; for k = £, r,, £ respectively. The appear- 
ance of the term k t is to account for the divergence terms in equations (2.11), 
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(2.14), and (2.15) in which the vector w explicitly appears. This can be further 

explained by considering the continuity equation, the treatment being same for 
the other equations. The continuity equation can be rewritten as: 

^ + div [ q(u 4 w) ] =0 (2.35) 

which in curvilinear coordinates becomes 

^£l + JL[ Jg( u 4 w) ' a k ] =0 (2.36) 

dt dx k ~ 

It can be verified that: 

k t = w • a k (2.37) 

and 

K' = (u + w) ■ a k (2.38) 


In the formulation of the continuity equation in the absolute frame the term k t 
appears in a strikingly similar manner and denotes the contravariant compo- 
nent of the grid velocity vector. In the present formulation, however, there is no 
grid motion involved and k t should not be associated with any grid velocity. It 
is just the contravariant component of the vector, w = — Q x r. The expres- 
sions for transformed shear stress and heat flux terms are as follows: 


T kx — k XX "t kyTyx 4" k z tzx 

T ky = k x txy 4 kyTyy 4 k Z T Z y 

T fa = k x TxZ 4 kyty z 4 k z T zz 

Qk = uT kx + vT k y + wT kz + + k yQy + 


(2.39) 


where 
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T rr ~ 


--U 

Re 3^ 


2(^f + ^ + + ^ + ^ 


du A _ y du, (t dv J. „ dv x >. aU) 


/ t dW , „ dw , * dlVs 

(?2 s| + + ?2 af* 


T >y ~ 


——u 
Re 3” 


Ort J. n 0V -U f OU ) (t OU , OIZ , t £U) 

2( + Vy dij + * ap ~ (£ * a* + TJx dr, + %> 


dv 

drj 


dtA 

a£ J 


du 

d£ 


du 


du. 


(t dw + „§!£ 4. f dW) 
V 5 z it + * 1 * An + a?"' 


3| 


dt] 


% 


T 7.7. 


±Zu 

Re 3” 


9/t 9w , dw , t dw, _(tdu + du + tdu ) 

mz li + r]z d V + a£ a£ a»/ e * ar 

('t ay , „ ay , a^N 

( ^af + + 


T xy ~ T yx Yie^ 


It au , „ dM , f- 0 M\ I /-fc . 01/ . t 

( £y^t + ^^7 + ^1F + ( « x ^ + ^a» + 


a« 


a« 


dv 


dv 


dv, 


d£ 


drj *y% 


a| 


dt) * x d£ 


lyz ~ T zy ~ p^A* 


✓t ay , „ ay , y dv^ , /fc aw , dw_ , *. aw , 

<** « + ’*a* + 5z af ) + ^ af + ^ + a? ' 


Tzx ~ Txz " R e ^ 


/'t aiy , „ dw , v aty, , r t du , du . t du 
(SxTT + »/x^r + Cx^r) + + >7z3r + 


a£ 


a>7 


a?' 


a£ 


<7* = 
<7v = 

Qz = 


J. 

Re 

f* 

(y - l)Pr 

(fc £i_ 

( * x a£ 

1 

M 

it dT 

^ a£ 

Re 

(y - l)Pr 

1 

M 

(fc dT 
^ 2 d£ 

Re 

(y - l)Pr 


dT+t dT) 

c drj ^ x d£ J 

dT j. t dT\ 
'dri + Qy d r 

dT,fdJ\ 
drj + ^ z d£ ; 


(2.40) 


(2.41) 
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2.6.1 Thin-Layer Approximation 

Since the numerical simulation of the full Navier-Stokes equations are 
very computationally intensive, the thin-layer approximation is used in this 
study. This approach is useful because it can handle flows which have both invis- 
cid and viscous regions and interactions between the two [2]. The thin-layer 
approximation makes use of the fact that most grids for Navier-Stokes equa- 
tions are packed in the direction normal to the walls, whereas along the stream- 
wise direction they are usually very coarse and good enough for inviscid calcula- 
tions only [2]. In this study, § is the streamwise, rj the spanwise, and £ the 
pitchwise directions. Thus, all the viscous and heat flux derivatives in the 
streamwise direction are dropped. In the other two directions, only the normal 
second derivatives are considered while the mixed derivatives are neglected. 
This is done because they are insignificant in comers and lead to expensive cal- 
culations [2]. Applying the above approximation to equation (2.30) gives the fol- 
lowing final form of the governing equations: 

dQ . dF , dG , dH _ dGt, dHy o (2 42) 

dr d^ dTj d£ drj d£ 

where 

' Q ‘ 

QU 

Q = J Q v (2.43) 

gw 

6t . 

' gU' ] gV' ] T gW' 

guU' + guV + rjxp guW' + £*p 

F = J QvU' + £ y p , G = J Q vV ' + VyP , H = J Q vW + ZyP (2.44) 

gwlT + qwV' + rjzp gwW' + £zP 

e t U'+ P U \e t U'+pV J [e t W’ + pW 
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(2.45) 


(2.46) 


The terms U\ V' , W', and U, V, W are defined in equation (2.32). The ex- 
pressions for the transformed shear stress and the transformed heat flux terms 
appearing in the viscous flux vectors are: 


T kx = k x T xx + kyfyx + kz^zx 

Tfoy = k x T xy + kyTyy + k Z T Z y 

T = k x TxZ "t kyTy z 4" k z X zz 

Qk = uTfa + vT 'ky + wT + kxfjx + kyQy 4- k z q z 


(2.47) 


In the equation (2.47) above, k = tj, £ for flux vectors G v , H v respectively. 
Applying the thin-layer approximation in the k direction the viscous stress and 
the conductive heat flux terms will become: 
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r xy - r yx ~ 


Ty Z — T Z y 


Te» 


dU . l dV 

'dk X dk_ 


k dv + y dw 
* z dk ky dk 


r zx - r xz - Re v 


dw 


du 


h ^ + h 

* X dk + * Z dk 




( k — ) 
Kkx dk’ 

qx Re 

(y - l)Pr 

„ _ 1 


(y dk’ 

Qy Re 

(y - l)Pr 

„ -A. 


(k — ) 
{kz dk’ 

Qz Re 

(y - l)Pr 


( 2 . 49 ) 


( 2 . 50 ) 


The approach described in this chapter for the development of the govern- 
ing equations was not initially employed for this study. The initial route taken 
was via the cylindrical coordinates. However, it was later realized that the ap- 
proach described above is much more elegant and less time consuming. It also 
entails the minimum amount of changes to the existing code for it to work for 
this particular formulation. 

For completeness, the development of the governing equations using cy- 
lindrical coordinates in a rotating frame is presented in Appendix C. The final 
form of the equations derived in Appendix C are then transformed from cylindri- 
cal to Cartesian coordinates to see if they match the governing equations derived 
here in this chapter. 



CHAPTER III 


NUMERICAL SOLUTION METHOD 


The finite volume approach is utilized to obtain the numerical solution 
of the thin-layer Navier-Stokes equations which were derived in the previous 
chapter and are presented in their final form in equation (2.42). 



The curvilinear coordinate transformation maps the entire physical do- 
main to a computational domain which is discretized into collection of cells of 
unit volume. Thus, A£ = Arj = At, = 1. The finite volume formulation is ob- 
tained by integrating equation (2.42) over a unit volume computational cell, 
which is shown in the figure below. 
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In the cell-centered finite volume approach the dependent variable vector Q, is 
assumed to be the flow quantities at the cell center, and is uniform throughout 
the cell. The fluxes F, G, H, G v , H v are assumed to be uniform over each cor- 
responding cell face. The integral form of equation (2.42) after spatial discreti- 
zation, can then be written as: 

+ KG. i - G. ,) - (G v . - G v ,)] A&C (3.1) 

2 J 2 . 7*2 J 2 

+ [(H k+ , - H k J - (Hl^ - H v k J]A^AT} = AZArjAt; S 

Defining the central difference opera tor das <3, ( ) = ( ). i - ( h^andnot- 

1 + 2 1 2 

ing that the finite volume cell has sides of unit length, equation (3.1) can be re- 
written as: 

^ + d:(F) + 6 AG - G v ) + d k (H -H v ) = S (3.2) 

OT 1 J K 

In the above definition of d, l + b represents the face adjacent to the cell center 

in the positive direction of l and l ~ b represents the face adjacent to the cell 

z 

center in the negative direction of l, where l = i, j, k for£, rj, £ directions re- 
spectively. 


3.2 Implicit Formulation 

An implicit scheme is utilized to numerically integrate equation (3.2) be- 
cause these schemes can handle large CFL numbers. In problems involving vis- 
cous computations, very fine grid spacing is chosen to resolve the viscous effects, 
due to which the tolerable minimum time step cannot be too small for the code 
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to have any practical value [2] . A first-order time-accurate implicit scheme to 
numerically solve equation (3.2) can be implemented by the following algebraic 
equation: 

AQ n = - Ax 0 R n + 1 ) (3.3) 

where 

A Q n = Q n + 1 - Q n 

R n + 1 = d- (F) n + 1 + dj (G) n + 1 + d k ( H) n + 1 - (3.4) 

dj (GT +1 - d k (H v ) n + 1 - S n + 1 

As can be observed from the above equation, the convective and diffusive fluxes 
and the source term are all treated implicitly. 


3.3 Newton’s Formulation 

The Newton’s iteration method is used in this work to solve equation 
(3.4). A detailed explanation of the method is presented in [2] and [7] . For a sys- 
tem of nonlinear equations the Newton’s method can be stated as [7]: 

F'(x k ) (x k + 1 - x k ) = - F{x k ) (3-5) 

Dividing equation (3.3) by J, will result in 

Aq n = - Ax CR” + 1 ) (3.f>) 


where 



The Newton’s method can now be applied to the following vector equation: 

L(q n + 1 ) = 0 (3-3) 


where 
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L(q n + 1 ) = q n + 1 - q n + At (R n + 1 ) (3.9) 

The resulting matrix equation is: 

[ /-Jf +Z¥(M )*’ k ~ 1 ] Aq k ~ l = 

dq (3.10) 

— [ q k ~ l - q n + At (f?*’* _1 ) ] 


The operator ■ indicates that the difference operator M acts on the delta form 

of the unknown vector Aq k ~ l , where Aq k 1 = q k — q k 1 [2]. It can be ob- 
served from equation (3.10), that the time derivative term is present in the resid- 
ual which arises by applying the Newton’s method to the nonlinear system of 
equations. The sub-iterations help in obtaining a better approximation of the 
solution at each time step. This method, is therefore, better for unsteady prob- 
lems. The various operators and terms appearing in equation (3.10) are defined 
as follows: 


• = <5 f (A)* ,a_1 ■ + dj ■ + d k ( C )*-*- 1 • 

dj (. B v )*’ k ~ 1 • - 6 k ( C ")*’*- 1 • 

i a dF p dG dd dG f'V 

where A = B = C - B - , C 


R *,k-i = 6i + dj (G)*’* _1 + d k (H) 

h 


*,k-i _ 


6 (G v )*’ k ~ 1 - a. (. H v )*’ k ~ l ~ S**- 1 

J 


dH v 

^ (3.11) 


The superscript * is used in the above equation to specify the time level 
for computing the metric terms. Janus [13] found that in order to maintain the 
mass flow without discontinuities when crossing the interface between a rotor 
and a stator, the metric terms at old time level n and q° should be used for the 
first iteration. For the subsequent iterations, the metric terms at the present 
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level n + 1 and q k ~ 1 are used. A, B, and C are the flux Jacobians, which are 5 
x 5 matrices with real eigenvalues. The form of the flux vectors derived in chap- 
ter 2 is identical to the flux vectors obtained in the absolute frame formulation. 
The derivation of the transformed flux Jacobians and their eigensystem is given 
in detail in Janus’s thesis [8]. 


3.3.1 Flux Vector Splitting 

The flux Jacobians have real eigenvalues which correspond to wave prop- 
agating speeds in each curvilinear direction. The sign of the eigenvalues dic- 
tates what information should be used for the evaluation of the fluxes. To incor- 
porate this approach the flux vectors are split using the Steger- Warming flux 
vector splitting method described in [9]. The flux vectors are split as follows: 

F = F + + F~ 

G = G + + G~ (3.12) 

H = H + + H~ 


where the superscript of plus denotes the flux sub-vector corresponding to the 
positive eigenvalue so that this flux sub-vector is computed using information 
from the negative direction. Similarly, the flux sub-vector with the superscript 
of minus is evaluated using information from the positive direction. Utilizing 


this concept a new set of flux Jacobians can be formulated as follows: 


dF + 

B + = 4^ + > 

SI' 

II 

+ 

dq » 

dq ’ 

dq 

dF 

B~ = — , 

C- 

dq ’ 

dq ’ 

dq 


(3.13) 


The detailed derivation of these new flux Jacobians is available in Belk’s dis- 
sertation [10]. 
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3.3.2 Viscous Flux Jacobian 

The viscous flux Jacobians B v and C v , defined in equation (3.11) are 
computed numerically in this study. The numerical derivative used to obtain the 
elements of the Jacobian matrix is [11]: 

F[ (x + he m ) — ^ 

a im ~ h 

where e m is the m tfl unit vector and h ~ JmachXne e. As will be explained lat- 
er in this chapter, the viscous fluxes at a face are functions of the right and left 
dependent variable vector. Thus, the viscous numerical derivative will consist 
of two parts which can be written separately as: 


&v + 

= £L 

(q j + he m , Qj+\) Q.j + 1 ) 


h 

/ N 

to 

i 

= 2L 

(qj, qj + 1 he m ) — G v t (q qj + 1 ) 


h 

< C V + 

. n 



h 


.3 1 

Qk + i he m ) ~ Hf ( qk> Qk + \) 


h 


Equation (3.10) can be reformulated by using the new flux Jacobians as: 

[ / - 27 (l^)*’*' 1 + M {(M + )*’ k ~ l ■ + • /] Aq k ~ l = 

oq 

- [ q k ~ l - q n +Ai ( R ] 


where 


28 



= <3; (A + )*’*- 1 • 

+ dj (B + )* ,k ~ 1 • + d k 

(C + ) 


- dj [(B v ) + ]*’ k ~ l 

• - d k [(CV]*’*- 1 • 



= d t (A - )*’* -1 • 

+ dj (B~)*’ k ~ 1 ■ + d k 

(C-) 


- dj KB”)-]*’*- 1 

• - d k [(C v )-]*’ k ~ 1 • 



The difference expressions that have been mentioned till now are first order ac- 
curate in time. A general form of equation (3. 16) which incorporates a wide vari- 
ety of different time accurate schemes is written below [2] : 


1 + j _ 1 

6Az dq 


+ (M + 




Aq k ~ l = 


1 + \p 
OAr 


„*-i 


— ^(27 [ (0 - 1 )R n > n - 0 R*' k ~ l ] (3.18) 
1 + xp l 

-I- xp Zlq" -1 } 


A few important time accurate schemes are Euler implicit (0 — 1, xp — 0), 
three-point backward (0 = 1,^= V 2 )> and trapezoidal (0 = 1/2, xp = 0). 
The first scheme is first order accurate whereas the other two are second order 
accurate in time. The three point backward scheme has been used for the work 
carried out in this study. 


3.4 Flux Formulation 

The computation of the RHS of equation (3.18) requires a method for the 
computation of the fluxes. The term R defined in equation (3.4), consists of con- 
vective and diffusive flux terms, each of which is computed differently. 

3.4.1 Convective Flux Computation 

The convective fluxes are hyperbolic in nature which implies that the ei- 
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genvalues associated with the flux Jacobians are real. Thus, there is a preferred 
direction of information propagation in such systems. Thus, to be consistent 
with the physics of the problem one needs to use the method of upwinding to 
evaluate the convective fluxes. One such method is Roe’s flux formulation which 
is a flux difference split scheme [12] . The method entails solving an approximate 
Riemann problem exactly at each cell face. However, Roe’s method is first order 
accurate in space which is too dissipative for simulation of realistic problems. 
Higher order accuracy is obtained by adding corrective flux terms with limiters 
[2] . The limiters act as a switching mechanism to control the addition of the cor- 
rective flux terms such that the formulation is higher order in smooth regions. 
In regions containing shocks or contact discontinuities the formulation is first 
order accurate because in such regions higher order schemes are dispersive. The 
limiters used in the corrective flux terms are Roe’s Superebee, van Leer’s and 
minmod limiters. The detailed development of the above formulation is pres- 
ented in [13], [14], and [15]. 

3.4.2 Diffusive Flux Computation 

The evaluation of the diffusive fluxes defined in (2.45) requires the nu- 
merical representation of derivatives needed for the stress calculations which 
can be noted in equations (2.47) and (2.48). The metric terms such as £ y , ...., 

are also needed instead of J£ x , J % y , ...., which are the projections of the face 
areas in the Cartesian directions and are computed in the code. The diffusive 
fluxes are parabolic in nature so that there is no preferred direction of informa- 
tion propagation for systems. Thus, the central differencing scheme is the most 
popular and easiest way to compute the derivatives appearing in the diffusive 
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fluxes [2] . The implementation of the above approach is illustrated by consider- 
ing an example stress term, r**: 


Txx Re 3^ 


du 


dv 


dw 


o b — b i-ii _ k 

Zkx dk y dk 2 dk 


The following numerical formulae are used to calculate the various terms that 
appear in the above equation: 


( k x ) 


< Jk * >* +i 
¥ J k + J k +i > 


’ ( k y \+y ~ i 


( Jk y 


2 2 ^ + ^k + l ) 


< — ) 

( dk '*+t 


1 - U k + 1 U k , 


( — ) 

' 7 . 'I 


'*+± u *+i 


- 


(3.19) 


(3.20) 


The central difference scheme in equation (3.20) is of second order accuracy. 


3.4.3 Turbulence Modeling 

Equation (2.42), which is the thin-layer Navier-Stokes equation, becomes 
Reynolds-averaged thin-layer Navier-Stokes equation which is the time aver- 
aged Navier-Stokes equation for turbulent flows, when fi is written as: 

fi = n t + fi t (3 - 21) 

Similarly the thermal conductivity is written as: 

k = k t + k t (3.22) 


or 

k = t = Hi + tl 

Pr (y - 1) Pr z (y - 1) Pr z (y - 1) 
where Pr z = 0.72 and Pr z = 0.90 for air [2]. 


(3.23) 


The Baldwin and Lomax turbulence model [16] is used in this study for 
the calculation of n t . A detailed description of the model is given in [2] . This mod- 
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el is a two-layer algebraic eddy viscosity model in which fi t is evaluated differ- 
ently for the inner and outer layers. Since there is no clear demarcation between 
the inner and outer layers, / i t at each cell face is calculated by both methods and 
then the smaller value is chosen [2]. 


3.5 The N-Pass Scheme 

Equation (3.18) can be written in the general form as: 

( L + D + U)X = b (3-24) 


where L is strictly lower block triangular, D is block diagonal and U is strictly 
upper block triangular. The system of equations defined above can be solved us- 
ing the “N-Pass” algorithm. The “N-Pass” algorithm can be viewed as a relax- 
ation scheme based on the symmetric Gauss-Siedel algorithm [17]: 


( L+D)X 1 + UX°=b 
LX 1 + ( D + U ) X 2 = b 


1 st pass 


( L + D )X 2N ~ l + U X™- 2 = b 
L X™~ 1 + ( D + U ) X™ = b 


> N^ 1 pass 


(3.25) 


The nature of L, D, and U can be exploited to reduce the number of matrix vector 
multiplications needed to be carried out [7]. The modified two-pass scheme can 
be shown to be equivalent to one pass of the “N-Pass” scheme described above. 
This has been illustrated by considering a representative one-dimensional sten- 
cil in [2] . Since, more than two iterations are needed for the Gauss-Siedel itera- 
tive method to yield a converged solution, the “N-Pass” scheme allows the scope 
for a more accurate and stable solver. This is so, because with each additional 
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pass one can add two additional Gauss-Siedel iterations, which gives a better 
solution. In this study, only one pass has been used to obtain the solution. This 
scheme has also been vectorized using the diagonal plane processing developed 
in [13] and [10]. 

3.6 Boundary Conditions 

The characteristic variable boundary conditions as applied to the non- 
conservative vector form of the transformed Euler equations with source terms 
developed in [18] is used in this study. For the viscous computations, at the 
walls, the adiabatic wall, no-slip boundary condition developed in [14] is ap- 
plied. The inflow boundary condition that is developed in [2] which specifies the 
reservoir conditions for internal flows, called the total-condition- preserved 
boundary condition does not change between the two formulations. This is so be- 
cause inflow is assumed to be uniform at the grid entrance and only in the x- 
direction which is normal to the plane of rotation. The only difference in the 
boundary conditions that is made when one changes from the fixed frame ap- 
proach to the rotating frame approach arises in the subsonic outflow boundary 
conditions. The subsonic outflow boundary conditions with source terms are 
coupled with the simple radial equilibrium equation described in [19]. The sim- 
ple radial equilibrium accounts for the swirl produced downstream of the blades 
[ 2 ], 

3.6.1 Subsonic Outflow 


Subsonic outflow condition is illustrated by the figure below. 
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Figure 3.2 Codirectional Outflow 

In the above situation, depicted by Figure 3.2 four pieces of information 
come from inside the computational domain and one from outside. If the outflow 
boundary is a constant £ plane and the outflow is in the positive £ direction, then 
the characteristic variables that are evaluated from inside the computational 
domain are [18]: 


%X (6 - 4 ) + £* v 

c o 


~ §y W 


(3.26) 


(Q ~ + & v 

c o 


~ %y W 


+ r , 


Jin 


h - p> 

0 


- £z U + W 


Jb 


(y 

0 


(3.27) 


- Iz u + w 


+ In 


£z (Q - 4 ) + ty K - £x V 


Jb 


(e ~ 4) + £v u ~ £* v 


+ A 


<in 


(3.28) 
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pjVjj 

<?o c o 


+ £x U + £v v + £z w 


£o c o 


+ £ x U + £y V + £z w + r 4 


(3.29) 


where r = At - ^ * r , r = Pta M 1 S. The subscript “0” denotes that 

the matrix P is considered to be locally constant. The 1 and Af 1 are 
derived in [8] and S is defined in (2.46) so that one can write the elements of T 
as: 


r 4 = At [— £ z w Q — v Q] 
r 2 = At [£ x v Q] 
r 3 = At [£ x w Q] 
r 4 = At [- ty w Q + | 2 V Q] 


(3.30) 


By constructing the grid such that the outflow boundary is normal to the x axis, 
the metrics g y and £ z will be zero [19]. Utilizing this result in equations 
(3.26M3.30), the characteristic variables simplify to: 


P 

r, 

c 


p-4 


<3.31) 


'b - 
"b = 

f P 

{Q o c o 


oJb l c °k 


v in - zlr Wav & 

(3.32) 

w in + At V av & 

(3.33) 

+ “)„ = («£. + “L 

(3.34) 


where g 0 and c 0 are reference quantities. v av and w av are calculated using 
phantom points which are explained later in this section. The calculation of 
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pressure at the outflow boundary is carried out by using the result from the first 
three equations above in the radial equilibrium equation [19]: 


dp_ 

dr 


Q vl 


(3.35) 


A first order finite difference scheme is used to evaluate the derivative and the 
right hand side is calculated as the average of the quantity at the two points be- 
ing considered. This procedure will yield the following algebraic equation: 


P2 ~ Pi _ 1 

r 2 - r l 2 


te + a ‘) 17 


2 

01 


fe + <* 2 ) 

v ^ 
02 

\ 02 / 

r 2 



(3.36) 


where a = 


• -5 

o. 


and the subscripts 1 and 2 denote the points along the 


a vi 


spanwise direction. Letting f = 2 P and g — — one gets the following ex- 


r* T 

c o r 


pression: 


_ Pi + 1 ^(P 1 A+gl+g2) 
P2 " 1 - 


(3.37) 


The calculation proceeds from the casing (maximum value of r), where the pres- 
sure is taken to be the specified back pressure, to the hub. Once the pressure at 
the boundary is determined, the density and u can be calculated from the equa- 
tions (3.31) and (3.34) respectively. 

To implement the boundary conditions a layer of phantom points just out- 
side the computational domain is used. The conditions at these points are deter- 
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mined by the following relations: 

V'p = 2 xp h - v> in 


(3.38) 


where xp is any of the primitive variables and subscript p denotes a phantom 
point. Once the primitive variables are determined, the conservative variables 
can be calculated at the phantom points. Using the above definition, 
v av and w av appearing in equations (3.32) and (3.33) can be written as: 
v ; „ + v„ 

(3.39) 


Van 


Wav = 


w in + 


It should be noted that u p and w p in equation (3.39) are calculated using the 
solution at the previous time step. 



CHAPTER IV 


RESULTS 

The aim of the work carried out in this study is to show that, for turboma- 
chinery applications, solving the Navier-Stokes equation cast in the rotating 
frame has a significant advantage over solving the same in fixed coordinates. 
However, one has be sure that the change in formulation does not introduce any 
additional inaccuracy which will degrade the solution. The test case chosen in 
this study to validate the formulation is NASA’s Rotor 67, which has one row of 
rotating blades. 

All the computations were performed on the SGI Challenge 10000 XL 
computer. The calculations have been carried out with third order spatial and 
second order temporal accuracy. The higher order computations and viscous 
computations are delayed for a few cycles. 

4.1 Rotor 67 

The first stage rotor of a NASA Lewis designed two-stage fan, Rotor 67, 
has been used to validate the code for a single-rotating geometry in internal flow. 
The test case is a low-aspect-ratio, damperless, transonic axial-flow fan rotor. 
It has 22 blades of multiple-circular-arc design and no inlet guide vanes. The 
inlet relative Mach number at the rotor tip was 1.38 at the designed tip speed 
of 428.9 m/sec (16,043 rpm). The designed pressure ratio for the rotor is 1.63 at 
a mass flow rate of 33.25 Kg/sec [2]. 
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The grid used for the computations is very coarse. An H-type grid for the 
whole computational domain was generated by TIGER, an interactive turboma- 
chinery grid generation code developed by Shih [20]. The grid has dimensions 
of 55x31x25 with twenty one points on the blade axially and twenty six points 
in the spanwise direction. There is a gap clearance between the blade tip and 
the casing. Since this case is at zero angle of attack, only one blade passage 
needs to be simulated by taking advantage of the symmetry. 

Total-condition-preserved B.C. is used to apply constant uniform condi- 
tions at the upstream boundary. Axisymmetric B.C.s are used for re-entry 
boundaries. No-slip B.C.s are used on the blades, hub and casing. At the exit 
plane, the radial equilibrium boundary conditions with the source terms are ap- 
plied with the back pressure specified at the radial location corresponding to the 
casing. The Reynolds number based on the tip diameter at the designed speed 
and standard day conditions was about 8,000,000. 

The solution is considered converged when the following are true: 

1. The inlet mass flow rate changed less than 0.2% of the average of the max- 
imum and minimum values in 1000 cycles. 

2. Mass flow ratio (outlet / inlet) varied between 1 ± 0.005 in 1000 cycles. 

For the above test case three conditions were run with the non-dimen- 
sional back pressures of 0.75, 0.85, and 0.90 which correspond to choke, near 
peak efficiency and near stall conditions respectively. Since the problem can be 
viewed as a steady state flow in the rotating frame, the relative frame code is run 
with both, local time stepping and minimum time stepping. The absolute frame 
code can only be run with minimum time stepping. Figures 4.1, 4.3, and 4.5 show 
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a comparison of mass flow histories between the relative frame and the absolute 
frame codes, for the three conditions, respectively. The comparison of the behav- 
ior of the mass flow ratio is presented in Figures 4.7, 4.9, and 4.11. 

As can be observed from the mass flow histories, the absolute frame code 
requires a much longer time to yield a satisfactorily converged solution. The ab- 
solute code was run with a CFL=70 which corresponds to 2600 time steps per 
revolution, and with 2 Newton sub-iterations per time step. When using local 
time stepping the rotating frame code was run with a CFL=50; no Newton sub- 
iterations. While using minimum time stepping the relative code was run with 
a CFL=10000, corresponding to 18 time steps per revolution. However, the ini- 
tial 400 cycles were run with a CFL=500. Thereafter, the CFL was increased to 
10000. 

The absolute code requires 10400 cycles (4 revolutions) to produce a con- 
verged solution, whereas the rotating frame code requires a minimum of 1500 
cycles to give a converged solution when run using minimum time stepping. 
When local time stepping was used, a minimum of 1500 cycles were needed. 
Thus, there is no significant advantage of using the local time stepping method 
with the rotating frame code. This can be attributed to the very high CFL that 
can be used when running with the minimum time stepping method. It can be 
observed from the plots that as the back pressure is increased progressively, 
more number of cycles sire required to obtain a converged solution. The final in- 
let mass flow rates between the two formulations varies by less than ± 1%. The 
advantage can be better appreciated on a semi-log plot which has been shown 
in Figures 4.2, 4.4, and 4.6 for inlet mass flow and Figures 4.8, 4.10, and 4.12 
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for mass flow ratio. Another observation that can be made from the plots is that 
the mass flow ratio settles down much faster than the inlet mass flow for all the 
three cases. 

The speedup that one observes due to the above description is huge. The 
absolute frame code requires 38 hours of running time to complete one revolu- 
tion, whereas the relative frame code requires 10 hours running time for 2000 
cycles on an SGI Challenge 10000 XL. The CPU times required by the two for- 
mulations for the three different conditions are illustrated in the table below. 
Since there was no difference in the number of cycles required with the relative 
frame formulation when using local time stepping or minimum time stepping, 
the CPU times corresponding to the relative method using both approaches is 
the same. 

Table 4.1 

Performance Enhancement 


METHOD 

TIME 

(hrs.) 

back pr. = 0.75 

TIME 

(hrs.) 

back pr. = 0.85 

TIME 

(hrs.) 

back pr. = 0.90 

RELATIVE 

7.5 

12.5 

15 

ABSOLUTE 

152 

171 

190 


In order to ensure the accuracy of the solution, a comparison is made be- 
tween the relative Mach number contours obtained from by solving the two for- 
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mulations. The comparisons are made at two operating points, one near peak 
efficiency and the other near stall condition, and at 10%, 30%, and 70% span 
locations from the blade tip. The operating point near peak efficiency occurs at 
a non-dimensional back pressure of 0.85 which results in an inlet mass flow of 
33.89 kg/sec, which is 97.86% of the computed choked flow, which occurs at a 
non-dimensional back pressure of 0.75. For the same back pressure the absolute 
frame code has a mass flow rate which is 98.39% of the choked flow rate. Figures 
4.13-4.15 illustrate the comparisons for this operating point and for 10%, 30%, 
and 70% spanwise locations respectively. 

The near stall condition occurs at a non-dimensionless pressure of 0.90. 
The relative frame code yields an inlet mass flow rate of 32.49 kg/sec which is 
93.82% of the computed choked flow rate, for this condition. The absolute frame 
yields an inlet mass flow rate which is 94.91% of the choked flow rate. The rela- 
tive Mach number contours for this condition at the various spanwise locations 
are shown in Figures 4.16-4.18. 

It can be observed from the figures that there is a good agreement be- 
tween the relative Mach number contours. Thus, it can be concluded that the 
formulation presented in this study does not introduce any additional errors. 
The discrepancies between the computed relative Mach number contour plots 
and the experimental plots are the same for the two formulations and are not 
discussed in this work. 

A detailed description of the flow field and the agreement and discrepan- 
cies between computed and experimental results is given in [2]. Another ob- 
servation that was made during the course of this study is that the relative 
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frame code does produce a converged solution for a non-dimensional back pres- 
sure of 0.95 whereas the absolute frame code exhibits a stall condition beyond 
a non-dimensional back pressure 0.90. 




Inlet Mass Flow (K^sec) Inlet Mass Flow (Kg/sec) 
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Figure 4.5 Mass Flow History (back pr. = 0.90) 
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Figure 4.6 Mass Flow History (log-linear, back pr. = 0.90) 
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Figure 4.7 Mass Flow Ratio History (back pr. = 0.75) 
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Figure 4.8 Mass Flow Ratio History (log-linear, back pr. = 0.75) 
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Figure 4.9 Mass Flow Ratio History (back pr. = 0.85) 
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Figure 4.10 Mass Flow Ratio History (log-linear, back pr. = 0.85) 
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Figure 4.11 Mass Flow Ratio History (back pr. = 0.90) 
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Figure 4.12 Mass Flow Ratio History (log-linear, back pr. = 0.90) 
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CHAPTER V 


SUMMARY AND CONCLUSIONS 

The results presented in Chapter IV amply demonstrated that solving the 
Navier-Stokes equations in rotating coordinates is a worthwhile exercise and 
provides a huge advantage in terms of savings in computer time. The savings 
are due to two main advantages of the rotating frame approach over the inertial 
frame approach. Since the flow can be considered to be steady in the rotating 
frame, no Newton sub-iterations are required so that less computing time is re- 
quired for each time step. This is augmented by the improved stability of the ro- 
tating frame code which allows it to handle a far larger CFL than the absolute 
frame code. This reduces the total number of cycles that are needed to obtain a 
converged solution. Also the need to move the grid is obviated due to which the 
metric terms do not need recomputing after each time step. All the above factors 
contribute to the huge savings in computing time. The speedup in CPU time ob- 
tained for the test case in this study was a factor of 16. Moreover, the approach 
used in this study did not utilize multigrid method or Jacobian freezing. These 
two factors can contribute to further savings in computing time. 

The enhancement in the performance of the code is achieved by making 
a few minor modifications to the code to incorporate the rotating frame formula- 
tion. To summarize, they are the following: 

1. Freezing of the grid motion. 

2. Calculating the term k t , analytically as described in Chapter II. 
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3. Addition of the source term implicitly in the flux balance. 

4. Addition of the source term in the subsonic outflow CVBC which is 
coupled with the radial equilibrium equation. 

The results presented for the test case show that the solution obtained by 
the relative frame formulation is in very good agreement with the results ob- 
tained with the relative frame formulation. The significant improvement in sta- 
bility is useful for problems with different time scales in turbomachinery com- 
putation. Since the rotating frame formulation freezes the frequency associated 
with the grid rotation, the time scales associated with other frequency (e.g. blade 
fluttering ) can be used to investigate the unsteady behavior related to that par- 
ticular time scale. Thus, it can be concluded that this approach is the right direc- 
tion to take for simulation of flows through turbomachinery, which can be 
treated as steady state problems in the rotating frame. 

Future work will entail implementing the relative frame formulation for 
a multi-blade row configuration to simulate the flow through a stage. However, 
this would involve transformation of the dependant variable vector between the 
fixed and rotating frames before exchanging information between blade rows. 
Addition of the multigrid method and parallelization of the code are two other 
areas of future work. 
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DERIVATION OF THE VELOCITY OF THE MOVING COORDINATES 
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It was stated in Chapter II that the vector w can be proven to be negative 

of the velocity of the moving frame with respect to the stationary frame. The en- 
deavor in this exercise is to prove the above statement. Only the specific case 
encountered in this study is treated in this exercise. The figure below describes 
the situation by considering the plane of rotation. 



Figure A. 1 Plane of Rotation 

In the Figure A. 1 above, Y, Z represent the coordinate axes of the station- 
ary inertial frame whereas y, z represent the axes of the moving frame which is 
rotating with a constant angular velocity Q about the x axis, which points out 
of the paper; so that after time t the moving axes make an angle of Qt with the 
respective stationary axes. If one considers a fixed point in space denoted by the 
dot in the above figure, it will have fixed coordinates in the inertial frame, viz., 
(. X , Y, Z). However, it’s coordinates in the moving frame, viz., (x, y, z ) will change 
over time and can be expressed by the following relations. 

x = X 

y = Y cos(Qt) -I- ZsiniQt) (A.l) 

z = ZcosiQt) - ysin(^) 

Therefore, the partial time derivatives of the moving coordinates are given by: 


= 0 
dt 


dy 


— = Q [ ZcosiQt ) - YsinCf^)] = Qz 
dt 


= -Q [ycos(^) + Zs\n(Qt)\ = — Qy 
dt 

The vector w is defined in equation (2.3) as: 


dx l 

« dt 
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(A.2) 

(A.3) 

(A.4) 


(A. 5) 


According to the coordinate system described in Figure A.l, 
x 1 = x, x 2 = y, x 3 = z and for simplicity let i, j , and k represent the Car- 


tesian base vectors of the moving frame. Thus, 

w = i + ~zrj + = Qz j - Qyk 

dt __ dt J dt _ J J __ 


(A.6) 


The velocity of the moving frame with respect to the inertial frame is given by: 
i j k 

Q x r = | £) o 0 = — & Z J + (A.7) 

x y z 


It can be observed from equations (A.6) and (A.7) that: 


w = - Q x r 


(A.8) 



APPENDIX B 

CURVILINEAR COORDINATE TRANSFORMATION OF THE NAVIER- 

STOKES EQUATIONS 
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The governing equations in conservation law vector form using nondi- 
mensional variables in rotating Cartesian coordinates is given by: 


dq df dg dh = W dg^ , dW 

dt dx dy dz dx dy dz 


(B.l) 


where the vectors q, f, g, h, f, g v , h v , and s are all defined in equations 
(2.17M2.20). A general nonorthogonal, steady, curvilinear coordinate system is 
introduced in the ( x , y, z) space as follows: 


t = £(*,y,z) 

V = rj(x,y,z) 
t = t(x,y,z) 
r = t 


It should be noted that through equations (B.2) the rotating Cartesian coordi- 
nates are transformed to rotating curvilinear coordinates. In other words, the 
transformation takes place between two coordinate systems within the rotating 
frame. Using chain rule and writing the result in matrix form yields: 


r p i 


_ 

1 


■jr 

u 

dt 


1 

0 0 0 


dr 

d 

dx 


0 

T]x 


d 

d£ 

d 

— 

0 

Vy 


d 

dy 



dy 

N ^ 

1 


0 

Vz Cz 


d 


The inverse transformation in matrix form is: 
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■_a_‘ 

dr 


i 

0 0 0 


1 

d 

dZ 


0 

** ye z e 


d 

dx 

d 

— 

0 



d 

dT] 


x t) yt) z r) 


dy 

d 


0 



d 

Qj> 

r 

i 


X K ye \ 


dz 


This implies that: 


o 

o 

o 

i-H 


1 

►— 1 
O 
O 
© 

0 Vx Kx 


0 X S y ( Z i 

0 £y TJy 


0 Xfj y v z v 

0 £z Vz £ z 


I 

o 

H 

trr 

-rr 

N 

iTY 


(B.4) 


(B.5) 


Thus, J, which is the Jacobian of the inverse transformation, can be written as: 


J = 


10 0 0 
0 z % 

0 x rj yr) z rj 

0 y K z K 

l 


x e ye z s 

x tj yr) Zf) 

x t y$ z 5 


(B.6) 


The metric terms in equation (B.3) can now be expressed by the following rela- 
tions: 

£* = 'j(yi? z £ — z^) ~ £* = — 

= XfjZ^) t)y -j(XgZj- 2|JC^-) — -j(ZgXjj X^Zfj) (B.7) 

£z = »?z = j(y^ 5 - x^y ? ) £ z = - >$**) 
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The governing equations in curvilinear coordinates, therefore are: 

I + «.$«'-/■> + 

+ l,^(g - «*) + Vy£j<g - g v ) + Ky ^(g - g V ) 

+ i.4(A - A") + - A") + - A") 


= S 


Multiplying throughout by J yields the following: 


I +«.^v-n + fc^v-n + fc^<r-n 

+ 1,^(8 - g") + >/,&<« - *•) + z,&(e ~ g ° ) 


d>7 




+ - *") + Km 


-Hpih - h v ) 

oQ 


= Js 


Using the identity: 

- ^wj 

the equation (B.9) becomes: 

A(J 9 ) -q^ + jzU(,U+^ + ZM 1 + ^t%/ + w + iM <] 

+ 4 u<w + + w)i - </ - ^ wt«>] 

- (? - + gj(Jn, ) + ^w,)] 

- (A - A“)4(J|,) + Atfj.) + |(« 

= 4[J(&/' + ^ + l*A»)] + ^[J(^ + Vyg" + fcA')] 

+ 4 i>kw + iyg v + Kzh v )\ + js 


(B.8) 


(B.9) 


(B.10) 


(B.ll) 


Since the curvilinear transformation is steady: 
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= 0 (B.12) 

dr 

Using the expressions for the metric terms defined in equation (B.7) one can 
write: 


$<■*■> + + Ti {J ^ 

= - ^(x ( y,) + ^(x c y ( ) - -gjixpt) + - |(V{) 

= dyjfix_ , axj^y. _ jy a 2 * _ sx^y_ . dy a 2 * , ax a 2 y (b. 13 ) 

at atat drj atat at ata?7 at a^at at a?; at 

_ d y d 2 x _ dx a 2 y ay g 2 x a x a 2 y _ dy d 2 x _ a 2 y 

at a^at at a^at a^ atat at ata^ at ata^ a^ atat 

= o 


Similarly 

+ $«»*> + #*■> ■ 0 

+ !«,) = 0 


(B.14) 


Thus, the governing equations in curvilinear coordinates are: 
^ + dF , dG , 3 H = aF" aG” MU + q 

ar at a?7 dt at at 

where 


(B.15) 


' e ' 

pn 

pu 

pit; 

e t 


Q — Jq — J 


(B.16) 
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K = J( kxf + kyg + k z h ) 
qK' 

quK' + k x p 
= J I QVK' + kyP 
qwK' + k z p 
et K' + pK 


where K = k x u + k y v + k z w + k t 
k t = — k y Qz + k z Qy 
K = k x u + k v v + k z w 


(B.17) 


it = J{ kj 0 + kyg V + kgk V ) 
0 


= J 


M kx 

T ky 

Tkz 

Qk 


(B.18) 


S = Js = J 


0 ' 
0 

— qwQ 
qvQ 

0 


(B.19) 


In the above equations (B.17) and (B.18), K = F, G, H; K’ = U’ , V', W'; 

K = U, V, W; It = F\ G v , H v ; for k = £, rj, C respectively. The trans- 
formed viscous stress and heat flux terms in equation (B.18) can be expressed 
in terms of the derivatives with respect to curvilinear coordinates by applying 
the chain rule to the derivatives appearing in the viscous stress and heat flux 
terms of the viscous flux vectors in equation (2.19). The transformed viscous 
stress and heat flux terms are defined in detail in equations (2.39M2.41). 


APPENDIX C 

CYLINDRICAL COORDINATE FORMULATION 
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The Navier-Stokes equations in cylindrical coordinates, in an absolute 
frame with absolute velocity components using non-dimensional variables can 
be written as: 


dq df dg dh 

dt dz dr dd 


*F 


+ Ml + + k 

dz dr dd 


(C.l) 


where 


<7 = 


‘ rg ■ 

rgu z 

rgu r 

r 2 gu e 

re, 


(C.2) 


f = 


■ rgu z ' 
r(gu 2 + p) 

rgUrUz 

) r 2 gu e u z 
rgHu z 


, 8 = 


rgu r ‘ 

rgu z u r 

r(gu 2 + p) 

r 2 gu$i r 

rgHu r 


h = 


gu e 

gu z u e 

guru e 

r(gu% + p) 
QHu 6 


(C.3) 



' 0 ■ 


' o ' 


‘ 0 ' 


rr zz 


rc zr 


T Z0 

II 

rcn 

, 8° = 

rr rr 

h° = 

T rO 

r2r ez 


rT dr 
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rq z 


rq > _ 


Vo 
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gu 2 +p- x m 
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(C.4) 


(C.5) 


70 


In the above equations u z , u r , and u e are the axial, radial, and tangential ab- 


solute velocity components respectively. e t — — + \q{u z + uf + u'q) and 

VIZ 


H = Also, 


2 V \ 9 du z _ du r _ 1 du 0 _ Ur 
Xzz 3 Re I dz dr r d 0 r 

2 P \ 9 du r _ l du e _ du z _ Ur 
Trr 3 Re I dr r d 6 dz r 

- _ dUz du r | Ur 

Tee 3Re[ r d6 dz dr r 



_ J4_ 1 du r dUo _ Uq 

ri 6r Re I r dd dr r 

_ a \ Bu e , 1 du 2 

T ez T z0 Re dz r d0 


1 dT 

q z = r zz u z + x zr u r + r z6 u e + _ 1)Pr Jlte 

i r n ]dT 

q T - TrzU z + TrrUr + T^Uq + _ l)p r J ~dF 

l| /x I j QT 

q 6 = tteMz + r 6r u r + TqqUq + ^ — _ 1)Fr 


(C.7) 


For flows in a rotating frame, equation (C.l) can be cast in the relative 
(rotating) frame by the transformation [3] : 

= « re i - at < c - 8 > 

where Q is the rotational speed (negative with 6). Introducing equation (C.8) 
into equation (C.l) yields: 
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dq df dg d(h + Qq) df dg° a/T , r ( r a) 

It + te + aF + dd + ~dF + dd + k yKj * ] 

The above equation is in cylindrical coordinates but in the rotating frame. How- 
ever, the velocity components in equation (C.9) are the components of the abso- 
lute velocity, the only difference being that the components are now expressed 
in the relative frame. If one introduces a transformation in the relative frame, 
from the cylindrical ( z , r, 6) coordinates to the Cartesian ( X , Y, Z) coordi- 
nates which is of the form: 

X = z 

Y = rcos0 (C.10) 

Z = r sin 6 


equation (2.16) can be derived by repeating the procedure discussed in Appen- 
dix A. Using chain rule and writing the result in matrix form yields: 


'a' 

dt 


'1 0 0 0 ' 


' d ' 

dt 

a 

dz 


0 10 0 


d 

dX 

a 

dr 

d 


0 0 Y r Z r 

o o Y 0 Z e _ 


d 

dY 

d 

d0 



dZ 


The inverse transformation in matrix form is: 




1 

I — 1 

o 

o 

o 

1 


'a' 

dt 

d 


0 10 0 


a 

dX 




dz 

d 


N 

> 

O 

o 


a 

dY 




dr 

d 

dZ 


tsj ' 

<£> 

<£> 

o 

o 


a 

dd 


(C.12) 


From equation (C.10) the following results can be derived: 
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Y r = cos 6, Y e 
r Y = cos 6, r z 


— rsin#, Z r = sin0, Z e 
sin0, d Y = — — , 0 Z 


= r cos 6 

cos 6 
r 


(C.13) 


Thus, the Jacobian of the inverse transformation is given by: 

J = ryd z - rfly = f (C.14) 

Also, the following relations hold between the components of the absolute veloc- 
ity vector in the two coordinate systems 

— XL — 1^2 

Uy = v - u r cosd - u e sind (C.15) 

u z = w = u r sin0 + u e cos8 


Although, in the present context, the above relations give the relation between 
the components of the absolute velocity, it must be mentioned that they are valid 
for the components of any vector in the two coordinate systems. Applying the 
chain rule to equation (C.9) and using the relations derived above, yields the fol- 
lowing form of the equations: 


^ + Ml + W 

dt dX dY dZ dX dY dZ 


(C.16) 


where 



(C.17) 
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f = 


QU 


QU 2 + p 


QUV 

, g = 

QUW 


u(e t + p) 



QV ~ qQz 
qvu — quQz 
QV 2 + P — qvQz 
qvw — qwQz 
v{e t + p) — eflz 


, h = 


qw + gQy 
qwu + guQy 
gwv + gvQy 
qw 2 + p + gwQy 
w(e t + p) + efiy 


(C.18) 


r = 


0 

^ XX 

T X y 

z 

uXxx + vx X y + wx xz + q x 


g v = 


0 

T yx 

r yy 

T yz 

UX yx l UTyy 1 WXyz "i” q y 


(C.19) 


h v = 


UT , 


VT 


0 

X zx 

X zy 
X Z z 

+ 


zy 


WTzz + q z 


S 


0 ' 
0 

— qwQ 
qvQ 

0 


(C.20) 


It can be observed from the above equations (C.17MC.20) that the equations in 
Cartesian coordinates derived via the cylindrical coordinates are identical to 
equation (2.16). 



